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1 Introduction 



In lattice QCD, masses, energies and vacuum-to-hadron matrix elements are extracted 
from the large time asymptotics of Euclidean two-point correlation functions. Con- 
vergence to ground state matrix elements and energies proceeds with a rate of order 
exp(— A^"^) t) where A^"^-* is the energy difference of the first excited state and the ground 
state in the hadron channel characterized by a set of quantum numbers A. 

Hadron-to-hadron matrix elements of the type {A\h^\B) require us to consider in 
addition three-point functions, which contain two time separations, 

C(3)(t2,ti) = (O(^)(t2+tl)/lw(tl)[O('')]*(0)) (1.1) 

and the corrections to asymptotic behavior are 0( exp(— A^"^) t2), exp(— A^'^^ ti) ). Hence 
one wants both ti and t2 to be large. On the other hand, at large times the statistical 
errors of the Monte Carlo estimates typically increase: the noise-to-signal ratio grows 
with a rate exp((5(^)t2 + S^^^h), where 6^'^ is a positive energy difference. It is therefore 
necessary to compromise between the two sources of error. This compromise represents 
an important limitation to the achievable overall precision (statistical and systematic). 

As an example, consider the nucleon axial coupling qa- For this application, A^"^) = 
A^^) due to isospin symmetry and d^^'^ = J*-^^ « "^nucleon — ^rn-j^ \A,5j^ Due to the 
symmetry, ti = t2 is optimal in this simple but relevant example and 

• ti = t2 ^ l/A^^) w 0.5 fm is required to keep systematic corrections due to 
excited states small 

• but statistical errors become too large beyond a time t = 0(1/(5*^"^^) = 0(lfm)j^ 



Numerical results have been shown for this particular example in the reviews [8j[9| and 
recently in Refs. [3 , 10 . 



As a remedy one may try to reduce either the statistical uncertainties or the con- 
tamination by excited states. A general idea for reducing statistical fluctuations is 
to integrate over part of the configuration space analytically or by a multilevel algo- 
rithm |11|-(13| . In the pure gauge theory, a reduction of the growth of statistical errors 



as a function of time has successfully been achieved by multilevel algorithms 1 1 , 13] as 
well as by symmetry constrained Monte Carlo [l4,15 , but it appears difficult to make 
further progress in this direction for QCD with dynamical fermions. More radically, 
in lower-dimensional models a complete rewriting of the path integral led to simula- 
tion methods where errors can be kept constant at large time in specific channels |T6] . 
Returning to more moderate gains, in the Heavy Quark Effective Theory (HQET) one 
is in a special situation because 6^^^ is power divergent and depends strongly on the 

^Another example of interest is the B*Bn coupling where A''*' = A'^' due to heavy quark spin 
symmetry and the energy difference S^^^ = 5'^' is discussed in p 7]. 

^For our numbers we consider pion masses above the physical one. Close to the physical point the 
situation is somewhat worse. 
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discretization. An optimization of the discretization of HQET yielded a much reduced 
S^"^^ [7]. Despite these advances, we do not have a true solution of the signal-to-noise 
problem either in QCD or in HQET. It is therefore important to efficiently exploit the 
information present in an available set of gauge fields. In particular, due to translational 
invariance it is possible to construct volume- averaged estimators which should have re- 
duced variance compared to those in which one of the fields has a fixed position. This 
requires the stochastic estimation of the all-to-all quark propagator [l7| - |20| rather than 
the traditional calculation of a point-to-all propagator, and given the aforementioned 
exponential growth of signal-to- noise in Euclidean time, it is essential to "dilute the 
noise sources" (notation of [l9]) such that each has support only on a single time-slice. 

It is then natural to try to reduce the systematic corrections due to excited states. 
As a first step, one improves the interpolating fields 0^^\ usually by smearing (see 
Sect. [2]). In this way one may reduce the pref-actor of exp(— A^"^) t2). 



But in Ref. 21 it has been pointed out that by considering interpolating fields 
and the Generalized Eigenvalue Problem (GEVP), one can construct a time-dependent 
effective GEVP-optimized interpolating field where the gap A^"^-* is enhanced to 

A(^) = - . (1.2) 

Besides raising the relevant gap for the ground state n = 1 by a considerable amount, 
excited states n > 1 can then also be reached in each channel! This was demonstrated 
to work very well for a decay constant in HQET, i.e. a matrix element of the type 
(/|^w|0). In this work we show that it is also very advantageous for non- vacuum matrix 
elements. 

Moreover, we present a new formula, which involves the 3-point matrix correlation 
function and the GEVP eigenpairs. There is a summation over the intermediate time ti 
with t = ti + t2 held fixed. We consider now the "symmetric case" when initial and final 
states are related by a symmetry transformation (e.g. for g^), since the general case 
without the symmetry is more complicated as we will explain in the following sections. 
When A^^H ^ 1, the corrections to the matrix element are reduced 

from Cexp(-A(^)t/2) (fixed ti = t2 = t/2 GEVP) 

to C"A(^)texp(-A(^)t) ("summed" GEVP) , 

with some coefficients C, C given by matrix elements which are usually unknown. 
The N = 1 case of the general formula reduces to the summed 3-point function that 
has previously been used in early investigations of the nucleon sigma term [1] and 
QA 22 . The improvement of the convergence rate to the ground state has recently 
been emphasized in [2|[3]. Let us leave aside the pref-actors C, C", about which lit- 
tle can be said on general grounds. The remaining time-dependent factors satisfy 
A^"^) texp(— A^"^^ t) < exp(— A*^"^) t/2) for all t. Furthermore, when one is in the asymp- 
totic regime A(^) t > 1 the gain becomes significant: the "summed" GEVP method 
requires approximately half the total time separation for the same size systematic cor- 
rections. 
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The derivation of the formula for the matrix element, as well as the associated 
correction terms, proceeds roughly as follows. We start from the GEVP expression for 
the energy levels, in a theory with degeneracy i^^^^ = and augment the theory 

by a "source" term e/iw in the Hamiltonian. The matrix element {A,n\hyi,\B,n) is then 
obtained as a derivative with respect to e of the effective (time dependent) energy level 
at e = 0. This idea is worked out in Sect. [3j In Sect. |4j we report tests of the method 
for a toy model and also in a quenched QCD/HQET calculation. There we revisit the 
GEVP for energy levels, using the overlaps computed in the quenched case to fix the 
parameters of the toy model and examining the convergence of energies in the model. 
As we will discuss in the conclusions, we expect our method to be advantageous in a 
number of applications. First we set up the notation and describe the standard GEVP 
method of Ref. [2lj. 

2 Matrix elements from Euclidean correlators 

In this section we define the problem more precisely and describe the "standard" solution 
as well as the one using the GEVP. 

We want to compute a matrix element of a local operator /i^(x), 

Mmn = {A,m\h„{0)\B,n) , (2.1) 

where m,n > 1 label the excitations in each channel. The quantum numbers A and B 
associated with exact symmetries of the (lattice) Hamiltonian including e.g. momentum 
or flavors charges remain implicit in Ai. For the lattice Hamiltonian derived from the 
transfer matrix, we have H\A,m) = Em^^lAjUi) . We take the finite (space-) volume 
normalization of states {A,m\A,m) = 1, which is easily related to the relativistic one. 
The matrix elements are computed from correlation functions 

Cf(t2,ti) = {Oi^\t2 + h)h^ih)[of\o)r) (2.2) 
Cif\t) = (Oj^^(t)[of)(0)]*), dfit) = {of\t)[of\0)]*) (2.3) 

where Oj (t) are interpolating fields localized on a time-slice t with j enumerating 
different fields. They carry the quantum numbers B in the usual way. We now turn to 
different ways of reanalyzing the correlation functions. 

2.1 Standard ratios 

We consider m = n = 1 for describing the "standard" method, since it is largely 
restricted to ground states. One defines a ratio 

nt2,h) = ,^(^);; ' (^),' {(Efit) - Efmti - 1.)/2) (2.4) 
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for fixed i,j with 

t = t2 + h, Efit) = -dtlog{clf\t)), Efit) = -dtlogid^f\t)). 
Our lattice derivative is defined as dtf{t) = ^[f{t + a) — f{t)]. When the sectors (A) 



and (B) are related by a symmetry of the theory, the exponential factor in eq. (2.4) is 
unity, as Ef{t) = Ef{t). 

Many variations of the ratio are possible, e.g. replacing E^{t) — )■ E'^{ti). The 
ratio has a quantum mechanical representation (based on the transfer matrix of the 
lattice theory 1^ 

7^(^2, ii) = Mil + ci exp(-(^5^^ - E[^^)t2) + C2 exp(-(Ef ^ - E[^^)ti) + ... . 

(2.5) 

These correction terms have already been mentioned in the introduction. Note that 
replacing oj"^^ — )• X^fc '^\f^^k^^ O^^^ — )• a^jf^O^f^ with a specific choice of fixed 
coefficients a does not change anything in this formula except for modifying the pref- 
actors ci,C2. Instead, in the following section we turn to the use of the GEVP in order 
to change the exponential rates of the correction terms. 

2.2 Summed ratios 

An improved asymptotic convergence is provided by the effective matrix element 

A^f— = _a^a^7e(t-ti,ti) = 7Wii + 0(tAe-*^), (2.6) 

A = Min(£;5^) - e[^^ , E^^ - e[^^) . (2.7) 



Eq. (2.7) can be seen by explicit summation over ti of the transfer matrix representation 

in, 
1 



degenerate case E^jf'' = En^^ it has been used long ago 
gence rate has recently been emphasized in Refs. [2|[3]. 
non-degenerate spectra was introduced. 

2.3 GEVP improvement 



of eq. (2.4) and it is the = 1 case of eq. (2.15) (taking the limit to — ^ t)- For the 

and its improved conver- 



22 



n Ref. [3] the generalization to 



We here summarize Ref. j21j and apply it to the present case. We assume that we have 
N linearly independent fields Oj, with couplings to the low lying states. The labels A, B 



are dropped where statements independent of the channel are made. The GEVP 23 

{[CVn]i = T>f=lCij [Vn]j), 

C{t) Vn{t, to) = Xn{t, to)C{to)Vn{t, to) , (2-8; 



^For simplicity, we everywhere neglect terms which decay exponentially with the time extent of the 



lattice. 
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constructed from the matrices C*-"^^ , C^^^ at times t > to yields effective energies 24 

Ef{t,to) = -dtlog{Xn{t,to)) (2.9) 

which converge as [2l] 

Ef{t, to) = En + 0(exp(-Ajv+i,n t)) , Aiv+i,„ = En+i - En (2.10) 

provided one takes to > t/2, which we use here|^ 

The starting point for computing matrix elements is an operator (in each channel) 
which satisfies |21| 

Af{t)\0) = |n) +0(exp(-Ajv+i,„t)). (2.11) 

With the definitions 

N 

Vnit) = Vn{t + a,t), {u,w) = '^U*Wi (2.12) 

i=l 

Rn{t) = {vn{t),C{t)vn{t)r^'^eicp{Ef{t + a,t)t/2), (2.13) 
the explicit construction of A':n{t) is given by 

[•Af{t)]^ = e-^*i?„(t)(^;„(t), at). (2.14) 



With respect to 21 we have here made a specific choice for the relation of to ^-nd t, 
denoting the resulting Vn as Vn{t) with a single argument. 
We can then obtain the desired matrix element 

>^rn„ = >tfjt2,ti) + 0(exp(-A(^j_,i„ti),exp(-A5^j+,„t2)) (2.15) 

from 

M±{t2M) = (0|[if(t2)](^)/iw[[if(tl)](^)]t|0) 

= {v^^\t2),C^^\h,tMn\h))R^^\t^)Ri^\h). (2.16) 

Here we have reintroduced the labels A, B. Eq. ( |2.16| ) reduces to eq. ( |2.4[ ) for Na = 1 = 
Nb-, but taking Na,Nb larger improves the convergence and enables access to excited 
states. 

As before, one can formulate a simpler effective matrix element when (A) and (B) 
are related by a symmetry and only the m = n matrix elements are required. The 
symmetry means 

£;f'«(t,to) = Ef'^{t,to), vi^Ht,to) = vifHt,to) (2.17) 
for all t, to and n. The ratio (remember that we use the shorthand (t) = (t+a, t)) 

satisfies eq. (2.15) as well but may have reduced statistical errors. The leading error is 
minimized by the choice t2 = ti. 



*For fixed tg one has E^{t,t(,) = _E„ + 0(exp(min(A„+i,„, A„,„_i)t)) instead 



24 
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3 Improved method: sGEVP 

Here we combine the improvement by summation of Sect. 2.2 with the GEVP of Sect. 



2.3 



3.1 Symmetric case 



We consider the symmetric case eq. (2.17) and drop the labels A and B. As derived in 
App. 1^ the effective matrix element 

Ki^{t) = a Cff Un = Vn{t,to) (3.2) 

converges to the exact matrix element as 

Mf/{t,to) = Mnn + 0{AN+l,nteM-^N+l,nt)). (3.3) 

The formula assumes to ^ t/2 and the exact size of the corrections does in general 
depend on how we choose to, e.g. to = t — a vs. to/t = fixed. We shall demonstrate in 
Sect.[4]that the corrections in eq. (3.1) are very small generically. The label "s" stands 



for summed, since Kij{t) is a 3-point function summed over one argument. 
3.2 Asymmetric case 



In the situation when eq. (2.17) is not satisfied or if W6 WcLiit ci matrix clGniGnt .A^r 



with n 7^ m, we first define the estimator for the difference — 

E(t,to) = Ef'^it,to)-Ef'''it,to)'^^^Ei^^-E^^^ (3.4) 
as well as the energy shifted correlation function 

A,(t,to) = e-*^(*'*")c(f)(t). (3.5) 
The summed three-point function is then defined by 

K,,it,to) = a^e-(*-*i)^(*A')4')(t-ti,ti). (3.6) 

Everywhere we take ^ An approximation to the matrix element is 

\{Umit,to), [K{t,to)[Xi^\t,to)]-^ - Kito,to)]Wn{t,to))\ \ 



M':^^'{t,to) = -d, 



with 



^ [{U^{t, to), D(^){to)Um{t, to))iWn{t, to), Ci^) {to)Wn{t, to))^/^ j ' 

(3.7) 

D{t)u^{t,to) = Xm{t,to)D{to)Urr,{t,to) , (3.8) 

C7(^)(t) «;„(*, to) = A(f)(t,to)C(^Hio)w'n(t,to). (3.9) 
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We have observed numerically that in the case of tg) 7^ 0, it converges as 

Mt^^'{t,to) = >^„^„ + 0(A^exp(-A^o)), (3.10) 
see Sect. |4j The gap A is given by the minimum one in the two channels, 

A = Min «Vi - ^^^^ ' <Vi - • ^^-''^ 

Since the exponential convergence is now governed by to, there is no obvious advantage 



compared to eq. (2.16) unless one takes to ^ t. If statistical precision is good enough 
to allow for such large to, 3-point functions and 2-point functions with a maximal 
time extent of t are sufficient to obtain a convergence rate of 0(A t exp(— A t)) as in the 



symmetric case. Eq. (3.10) has not been proved formally, but our numerical investigation 



of toy models leaves little doubt that it is correct. 
4 Demonstrations 

We carry out two sets of demonstrations of how the various estimators for matrix el- 
ements work. First we consider toy models, prescribing spectra and matrix elements 
and do not take statistical errors into account. We construct "difficult" (large correc- 
tions due to excited states) and "easy" toy models. The second set of experiments 
is a quenched computation of the i?*-B7r-coupling, where realistic statistical errors are 
present. 

4.1 Models 

4-1-1 Definition of the models 

We first specify the spectra in dimensionless form. Two different ones are used below, 

roE^^ = n , roE^y = 1.1 x n . (4.1) 

The length factor ro is in principle arbitrary, setting the overall scale of the theory, but 
we think of it as ro ~ 0.5fm. Level splittings of around 1/ro are realistic in QCD, as 
the particle data book and lattice computations show. 
Next the overlaps 

^Pin = {0\Oi\n) (4.2) 



need to be fixed. In our HQET applications (see Sect. 4.2), we use spatially smeared 
quark fields to construct the fields Oj. We computed their overlaps ipm for n = 1, . . . , 5 
and i = 1, ... 7 using the GEVP "creation operator" [^^^(i)]^. For details we refer to 
the following section. Here we just take the approximate matrix 

/0.92 0.03 

tp^ = 0.84 0.40 0.03 -0.06 0.00 | (4.3) 
\0.56 0.56 




7 



corresponding to smearing levels 1,4,7 which is typically done in practice [25]. We 
observed a strong decay of the overlaps ipf^ with increasing n which suggests that a 
truncation with ipin = for n > 5 is realistic at reasonable time separations of the 
correlation functions, say t > ro/2. In any case, what we discuss here remains a model, 
but we expect it to be quite realistic. 

The matrix tj;^ represents a relatively comfortable situation which we may not 
always have. For that reason we also construct a more challenging case 

/0.9 0.1 -0.1\ / -l/(3n2) \ 

^?n\n<3 = 0.8 0.4 0.2 V^^„|4<n<20= 2n-2-(2n)-3/2 (4.4) 

\0.6 0.6 0.5 / \ l/(n- 1) / 

with a slow decay in n. We set ■0m = for n > 20. 

With the model matrix elements (again we note that these are not completely 
unrealistic) 

Mnn = 0.7-^ , Mn,n+m = ^ for m > , (4.5) 

and assuming the sectors A,B to be related by a symmetry, 

C\f{t) = ^.^fM')>-^''' = Clf\t), (4.6) 

n 

the model is completely defined. In particular we have 

Cl^^(t2,t^) = ^Ce-^"'*^A^n.e-^™*M^^),*„.. (4.7) 

n,m 

We refer to this model as SlSl. Replacing t/j^ by V'^ defines the model ClCl and finally 
with ip^,En^ for channel A and ip'^,E'i^^ for channel B we define the model SlCh. In 
other words we have the following table. 



model 










SlSl 






pW 




ClCl 


0C 






p(0 

-C'n 


SlCh 


0S 









4-1.2 Energies from the GEVP 



The corrections of the effective energies extracted from the GEVP, eq. (2.9), compared 



to the exact energies is shown in Fig.[T] We see how tQ > t/2 accelerates the convergence. 
As expected CI is a more challenging situation with larger corrections. One also sees 
that at short time {t/ro < 2) the dependence on tQ is typically not very dramatic. This 
feature has been observed in a number of practical applications 25 , 26 . Still, it appears 



dangerous to rely on this in general. In the left hand plot, we also observe the difference 
of the GEVP and a standard effective mass (dashed-dotted line). Here C22 is shown 
(the corrections for Cn are quite a bit smaller). 
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Figure 1: Corrections of eq. (2.9). The ground state n = 1 is plotted in red and has the 
small corrections. The first excitation (blue) is above. Dotted lines are for to = ro/5, 
while dashed lines are to = t/2 and full lines = t — a. Shown on the left is the model 
SI and on the right CI. The dashed-dotted lines show the corrections of the standard 
effective mass of the correlator C22 which approaches the ground state energy. 

4- 1-3 Matrix elements 

Let us start with the easiest situation, the extraction of ground state matrix elements 
m = n = 1 in the symmetric case. These are shown for two models in Fig. [2j The 
labeling of the different estimates is as follows: 



"ratio" dotted, black eq. (2.4) with ti = t2 = t/2 



"summed" dashed, black eq. (2.6) 



"GEVP" dashed-dotted, red eq. (2.18) with ti = t2 = t/2 



"sGEVP" blue eq. (3.1) 



The scale of the y-axis covers a variation of 10%. On the x-axis in this and the following 
figures we have for each method considered the total time extent of the 3-point functions 
since in a MC computation this generically governs the statistical accuracy. The graphs 
illustrate that the improved asymptotics of the sGEVP estimate (compared to the GEVP 
and the single operators) {N = 1) go hand in hand with smaller corrections at moderate 
time separations, t ~ rg . . . 2roj^ Among the estimates which do not use a GEVP, the 
summed method is generically better, at least when t is not too small. 

Diagonal {m = n) matrix elements for the degenerate case are shown in Fig. [3] on 
the left. For n > 1 only GEVP and sGEVP can be used for a systematic computation. 
Even though the scale of the y-axis is enlarged, we observe that sGEVP also works 

® Recall that the gaps of the models are A„+i^„ = 1/ro. 
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rather well for determining excited states. Note that with a GEVP with = 3 states 
(as is used here), the convergence of the m = n = 3 matrix elements is rather slow, but 
we show them anyway for illustration. It is strongly recommended to use a larger in 
a real computation of if statistical errors allow. 



0.76 
0.75 
0.74 



0.71 
0.7 
0.69, 



1 


...sGEVP, ^=t/2 


\ k=2 \ 


sGEVP, ^=t-a 


\ 

\ 


- - GEVP 

0|^, summed 


\ 

\ 


0|^, ratio 






/', k=1 


1 1 2 3 4 5 


6 7 8 



0.76 
0.75 
0.74 



0.71 

0.7 
0.69, 





...sGEVP, ^=t/2 




sGEVP, ^=t-a 




- - GEVP 




0|^, summed 




0|^, ratio 













t/r„ 



Figure 2: Effective ground state matrix elements, model SlSl on the left and model 
ClCl on the right, both shown as a function of the total time-separation of the 3-point 
function. On the left side, the sGEVP estimates for to = and Iq = t — a can't be 
distinguished in the figure. For the non-GEVP cases we show two different interpolating 
fields, Ofe, A; = 1,2. 

On the right of Fig. |3| we show the matrix elements Mi2- Here the sGEVP means 



eq. (3.7) with the energy shifts. The improvement compared to the standard application 



of the GEVP, eq. (2.18), is present but is not as impressive as on the left side, where no 
energy shifts are needed. We do not show levels above n = 2 since there a larger GEVP 
would be recommended as we discussed for the diagonal case. 

Finally, consider the situation where the spectra of the A sector and the B sector 
are different, as in the model SlCh. Example applications are i? — )• vr transitions or 
elastic form factors with momentum transfer. On the left side of Fig. [4]the A^n matrix 
element is shown. We again observe an impressive advantage of the GEVP methods, in 



particular of sGEVP over the standard ratio eq. (2.4). On the right side of the figure we 



study Aii2, where eq. (2.4) is not applicable. In this particular case, the amplitudes of 
the corrections of the GEVP effective matrix elements are relatively small and interfere 
destructively. It therefore happens to be more accurate than sGEVP for a range of t. 

In conclusion, the study of the models shows that the asymptotic convergence 
formulae also provide a very good estimate of the relative advantages of the differ- 
ent methods at intermediate t. In particular, consider first the degenerate case. The 
comparison of the asymptotic behavior, A^nn''^(t, to) — -Mnn ~ 0(tA exp(— tA)) vs. 
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Figure 4: Effective ground state matrix element Ain (left) and A4i2 in model SlCh. 



-^rm(*0)^o) — -M-nn ~ 0(tA exp(— toA)) suggests that to ^ t is needed to reach similar 
accuracy in the two cases and indeed we find this generically to be the case. One then 
needs the 3-point functions at twice the total time separation in GEVP compared to 
sGEVP. In the non-degenerate case, the convergence is governed by to in both GEVP 
and sGEVP. Here a very significant improvement is the change from a standard ratio 
eq. (2.4) to GEVP or sGEVP, see the left of Fig. [4} The right side of that figure shows 
that sGEVP yields considerable further improvement over GEVP when a large to is 
chosen. 
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4.2 The vr-coupling in the quenched approximation 

In the static approximation for the b-quark, the i?*i?7r-coupUng is denoted hy g. It is 
a leading order low energy constant in the heavy meson chiral Lagrangian f27H29] . As 



such, it is of considerable interest for chiral extrapolations of lattice results, employing 
a systematic expansion in 1/rn-b and rn^ / {Sir'^ F^) . The bare matrix element i^ 

g = l{B'mAkmBl+{0)) , A^ix) = M^h^l^M^) (4-8) 



with \B^{0)) polarized along the fc-axis, see also |2||30]. Note that here we use the 
normalization of states {B{p)\B{p)) = (i?^(p)|i?^(p)) = 2L^, which corresponds to the 
non-relativistic one in the infinite volume limit. We do not include the renormalization 
factor of the axial current anywhere. 

Our interpolating fields for B and B* are related by the exact spin symmetry of 
the static approximation and are generated by gauge-covariant Gaussian wave functions 
inserted between the static and the light quark field. Such gauge invariant interpolating 
fields were introduced in Ref. [22 1. We use exactly the ones of Ref. |2j with width 
r^f/ro = 0.36, 0.51, 0.62, 0.71, 0.87, 1.01, 1.13 (eq. (2.5) of Ref. |2|). For the present 
demonstration we work in the quenched approximation and the light quark mass is set 



to the mass of the strange as in 25 . An ensemble of one hundred gauge configurations 
is used on a 32 X 16^ lattice with spacing a ~ O.lfm and statistical errors are kept 
small by an all-to-all method |19| in combination with the static action "HYP2" [7j 
as done previously [2]. Here we use one hundred fully time-diluted noise sources per 
configuration. 

4-2.1 Approximate overlaps 

We first pick five fields Oi from our set with r^{/ro = 0.36, 0.51, 0.62, 0.71, 1.13. With 



the operator eq. (2.14) we can then compute the overlaps 



i^in = {0\Oi\n) = i;in{t) + 0{exp{-{EN+i - En)t) (4.9) 

V'in(t) = {0\diAf{t)\0) = Y,Cij{t)[Vn]j(t)Rn{t) (4.10) 

j 

where n = 1, . . . , 5 labels the excitations. The normalization of the fields Oi is irrelevant 
for all applications, but in order to have the interpretation of an overlap, we choose the 
normalization such that (7^(0) = 1. In this case a value of one for tpf^ means that 
Oi\0) = \n) without corrections, i.e. 100% overlap. Furthermore, we fix the signs by 
the convention ipin > for the value i which maximizes |V'm| at fixed n. 

Figure [5] shows examples of ipinit). Even if these are not precision determinations 
of the overlaps, they show interesting features. The local field shown in the first row 
has considerable overlap with all states considered. It is a bad interpolating field for 

®We thank Fabio Bernardoni for discussions on the effective theory and a check of the normahzation 
of p. 
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Figure 5: Overlaps for the static B-meson interpolating fields. Column n refers to state 
n, row 1 to the local (time component of the) axial current. Rows 2-4 correspond to 
the three interpolating fields used in our GEVP computation of g with radii r^f/ro = 
0.36, 0.62, 1.13 respectively. 



ground state physics. However the other fields with reasonable radii display a rather 
strong decay of the overlaps with growing n, indicating that the smeared fields provide 
a good basis of interpolating fields which couple little to excited states. Indeed, this 
figure demonstrates that these wave functions considerably reduce the overlaps to high 
excited states. Conversely, this also means that high excited states are difficult to access 
with these fields. 
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0.80 
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0.60 ; 



t/r„ 



t/r„ 



Figure 6: The matrix element ^ as a function of t/vQ. Left: sGEVP estimate eq. (3.1) 
with to = t/2, right: GEVP estimate with to = i/2. The error band is our best estimate 
determined previously with very high statistics [2]. 



Reading off approximate plateau values, we extract the model i/j^ for Sect. 4.1 , This 



13 



g (t) 



0.80 

:) ; 

0.75; 

0.70 ; 

0.65; 

0.60 ; 



0.80 

g (t) 

0.75: 

0.70 ; 

0.65; 

0.60 ; 



t/r. 



t/rn 



0.80 

g (t) 

0.75; 



0.70; 
0.65; 

0.60 ; 



t/r„ 




Figure 7: Estimates of 5 as a function of t/rQ. Left: summed ratio estimate eq. (2.6), 



right: ratio eq. (2.4). On the top the interpolating field with the biggest overlap with 



the ground state is shown (r„f = 1.13ro). The bottom two figures are for r^f = 0.62 tq. 



0.80- 0.80- 

g (t,) - g (t; - 

0.75 7 0.75 7 



0.70- 0.70- 

0.65 7 0.65 7 

O.6O F O.6O F 

0.5 1.0./, 1.5 2.0 0.5 1.0./, 1.5 2.0 

Figure 8: Standard plateau plot for g. The ratio TZ{t — ti,ti) is considered as a function 
of ti/ro with t fixed at t = 2.14ro. The latter is the value used in the so far most complete 
determination [31] , while earlier t/rQ ~ 4 was used on a lattice with a = 0.2 fm [32]. 
Left: wave-function with r„f = 1.13 tq, right: r^f = 0.62 tq. 
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model yields a qualitative understanding of the corrections. We believe that computing 
overlaps as done here may also be very useful for understanding the systematic errors in 
present extractions of nucleon matrix elements, namely the question of the magnitude 
of excited state contamination. Given an approximate knowledge of the spectrum, this 
contamination can be roughly estimated when the overlaps are known. Indeed, let us 
apply our approximate knowledge of 'i/'in < 0.1 for n > 1, together with the plausible 
assumption that matrix elements Aimn are of roughly the same magnitude as Aiu- 
Then, at time separations t = tq and for the best wave function, excited states make 
rather small corrections of order 0.1e~^, i.e. of the order of a few per cent. Therefore 
the matrix element g is a rather easy test case and all methods should be successful. 



4-2.2 The matrix element g 

In this section we show numerical results for g, computed with the various methods 
introduced above. The GEVP estimates for the ground state, displayed in Fig. [6j 
exhibit no corrections exceeding 2% once t = tq has been reached. For smaller sGEVP 
has smaller corrections than GEVP. However, at large times the statistical errors are 
increasing faster for sGEVP. 

For this particular matrix element and for the best interpolating field, the correc- 
tions for the summed ratio (Fig. [7| top) are somewhat larger than those of the standard 
ratio and again the summed method suffers from larger statistical errors at large times. 
However, in the case of a less optimal interpolating field (bottom of Fig. [T]), the summed 
ratio exhibits its superiority. 

For comparison, we also show the frequently used analysis where t is kept fixed 
(here at a value used previously in determinations of g [3l] ) and one looks for a plateau 
as a function of ti. With our precision one can observe the lack of a plateau in Fig. |8] 
for r^f = 0.62 ro, but with errors at a 1% level a false "plateau" would be observed for 
i/4 < ti < 3ti/4. This demonstrates the danger inherent in this method. The left hand 
side of the figure shows that for g a plateau with the correct height is obtained for a 
larger smearing radius. 



4-2.3 Excited state matrix elements gnm 

For excited state matrix elements. Fig. [9j only the GEVP estimates are applicable. 
They appear to work quite well for the first excitation and also for the second excita- 
tion a reasonable estimate can be obtained. The sGEVP again seems superior, as the 
deviations from our estimated asymptotic values are smaller. Figure 10 demonstrates 
these same features for an off-diagonal matrix element. 



5 Conclusions 



In this paper we have introduced the GEVP method with summation, denoted sGEVP, 
and we have examined several alternative methods for computing hadron-to-hadron ma- 
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Figure 9: The diagonal matrix elements 5„„ as a function of i/ro . Left: sGEVP estimate 



eq. (3.1) with to = t/2, right: GEVP estimate with to = t/2. An iV = 3 GEVP is used. 



The matrix elements are seen to be ordered gn+i,n+i < 9n 
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Figure 10: The off-diagonal matrix element gi2 as a function of t/rQ. Left: sGEVP 
estimate eq. (3.1 ) with to = t — a, right: GEVP estimate with to = t/2 — a. An = 3 
GEVP is used. 



trix elements. They have rather different asymptotic corrections due to excited states: 



ratio 

summed ratio 
GEVP 
sGEVP 
sGEVP 



eq. (2.4) 



eq. (2.6) 



eq. (3.1) 



eq. (3.7) 



eq. (2.16) 



exp(-A2,it/2) 
tA2,i exp(-tA2,i) 
exp(-AAf+i,it/2) 
iAAr+1,1 exp(-tAAr+i,i] 

tA exp(-toA) 



(just ground state) 
(just ground state) 

(equal energy case) 
(general case) 



In the last case, A is given by eq. (3.11) and one will typically use to = t/2. The form 



of the leading correction term of sGEVP is derived in the appendix for the equal energy 
case, while for the general one we deduced it from the numerical investigation of toy 
models. The GEVP correction term is known from [21] and for "ratio" and "summed 
ratio" it follows directly from the transfer matrix representation. 

We investigated two toy models constructed to be quite representative for heavy- 
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light meson matrix elements. In these models, the asymptotic forms of the corrections 
have been found to be a good guideline for the behavior at intermediate values of t, of 
the order t = {2 — 3)rQ. In particular we found that generically sGEVP has the smallest 
systematic errors, followed by GEVP. As a rule of thumb, sGEVP requires half the time 
separation of GEVP for the same systematic accuracy. 

A Monte Carlo computation of the B*Btt coupling g confirms our findings in the 
models concerning the systematic errors. In addition it allows us to make statements 
about the statistical errors which have to be balanced with systematic ones due to ex- 
cited states. Statistical errors grow more quickly as a function of t for sGEVP compared 
to GEVP, but a comparison at roughly the same amount of excited state contamination 
corresponds to a factor two between the values of t. A comparison at roughly fixed 
systematic error is shown in Table [T] We observe a minor difference for the ground 
state in advantage for sGEVP and the ratio of errors grows up to a factor five in the 
error for M33 for the considered matrix element. 

m n Mtn{t,t) Mt^{t,t/2) Mt^{t,t-a) 

1 1 0.004 0.003 0.003 

2 1 0.010 0.013 0.009 

2 2 0.032 0.012 0.013 

3 3 0.063 0.012 0.012 



Table 1: Statistical errors of various Gstiinators for Qrnn — 

M mn 

for t ^ ro. 

In the comparison of the different methods, one also has to consider the numerical 
effort to compute the effective matrix elements. We assume that one wants to control 
the corrections by computing the t-dependence of the estimators. In the summed cases. 



eq. (2.6) and eq. (3.1), this can often be done with a fixed number of quark propagator 
computations yielding a result for all t, by computing sequential propagators. The 
computation of g is such a case. In fact, since we have used a full all-to-all computation 
with "time dilution" (in the notation of [19]), also the GEVP estimate is obtained at 
the same expense. In contrast, if one only uses translation invariance on a time slice 



( "time-slice-to-all" ) , and for example varies t, keeping ti = t2 = t/2 in eq. (2.4) or 



eq. (2.18), then the required number of propagator computations is proportional to the 
number of t-values considered. In this situation the sGEVP method has an additional 
advantage. 

Taking statistical and excited state errors as well as the effort into account, sGEVP 
seems to be the overall most accurate, safe and efficient method. Given the difficulty 
in evaluating relevant correlation functions at large time separations and assessing the 
systematic errors, it still appears advisable to compare the different approaches in most 
cases. 

In our opinion the sGEVP method (and maybe the GEVP method) should be 
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applied to nucleon matrix elements such as gA or moments of structure functions, where 
large time separations are difficult to reach (3||8}{l0j and it is non-trivial to estimate 
possible contamination by excited states. In order to appreciate the last point, recall 
that in the standard ratio method, the systematic error drops like exp(— A2,it/2). In 
order to see such a term, one has to change t to t' such that the error term changes 
appreciably, say by a factor of three. One then needs t' — t ~ 2/A2,i ~ 1 frcj^ The 
summed ratio reduces this requirement by a factor of about two and the GEVP methods 
by a larger factor. This gains security in the detection of possible systematic errors. 

Acknowledgements. We thank Hubert Simma and Fabio Bernardoni for useful discus- 
sions and Ulli Wolff for valuable comments on the manuscript. This work is supported 
by the Deutsche Forschungsgemeinschaft in the SFB/TR 09 and by the European com- 
munity through EU Contract No. MRTN-CT-2006-035482, "FLAVIAnet". We are 
grateful to NIC and to the Norddeutsche Rechnerverbund for allocating computing re- 
sources to our project. Some of the correlation function measurements were performed 
on the PAX cluster at DESY, Zeuthen. 

A Derivation of the sGEVP method 



Here we give a derivation of the formulae of Sect. 3.1 



A.l Linear perturbation of the original theory 

We are here interested in the matrix element A4nn assuming the degeneracy of sectors 
A and B via eq. (2.17) and /iw(x)^ = /i„(x)j^ To arrive at an expression for the matrix 



element, we augment the original theory with Hamiltonian H (defined through the 
transfer matrix) by adding a perturbation term with strength e, 

^(e) = H + eKiO). (A.l) 

The twofold degenerate levels with energy E^^^ = E^^ = En are then split to ^^^(e). 
From standard degenerate perturbation theory one has 

E^ie) = En±eMnn + 0{e^) (A.2) 

with eigenstates |±,n) = [\B,n) it \ A,n)]/^/2 and 



(A.3) 



''We here again assume a gap of around 400 MeV. Close to the chiral hmit lower energy states with 
a gap of 2m^ exist, but probably have small overlaps with the typically considered interpolating fields. 

*The operator /iw(x) typically comes from the expansion of the electroweak hamiltonian density in 
terms of 1/Mw, but other applications are possible. For example the field h^(x) representing /iw(x) in 
the path integral may be h^{x) = {x) + (x) , with A^{x) = tp{x)T^'yk^stp, with the raising and 
lowering Pauli matrices in SU(2) flavor space. In this case, the matrix elements sought are the B* Btt 
coupling g or the nucleon axial coupling qa- 
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A. 2 GEVP in the augmented theory 



The desired -E^(O) is efficiently computed with a GEVP method as follows. We combine 
the interpolating fields o\'^\o^^^ from Sect. [2] 



N. 



(A.4) 

Since A, B correspond to different sectors (e.g. different flavours) we have {A, m\0\ |0) = 
= {B,n\o\^\Q) . Expanding the path integral to first ordei 
diately that the combined 2N x 2N matrix correlation function 

Q,(t,e) = (O,(t)O](0)) 

has a simple block structure, 



in e one then sees imme- 



(A.5) 



(C(^){t) em , , 



(A.6) 



up to first order in e. The entries C^^^ = C^^^ were defined in eq. (2.3) and K in 
eq. 



The generalized eigenvalues A„, eq. (2.8), determine effective energies 
Ef{t,to,e) = -dt log(An(t,to,e))- 



(A.7) 



In the augmented theory, an extra argument e has been added to A„ for clarity. Eq. (2.10) 
describes the corrections by which En{t,to,e) differ from the exact energy levels. Dif- 
ferentiating that equation with respect to e yields 



A^„„ + 0(Ajv+i,niexp(-AAr+i,ni)) • (A. 



e=0 



It remains to give an explicit expression for E'^ {t, to) in terms of the correlation 
functions, which is equivalent to a solution of the GEVP to first order in e. The 2A^ x 2A^ 
GEVP equation, C{t,e)vn{t,tQ,e) = Xn(t,tQ,e)C{to,e)vn{t,to,e), separates into the two 
independent ones 



[C(^)(t) ± eK{t)]ut{t, to, e) = A±(t, to, e)[C(^Hio) ± eK{to)]ut{t, to, > 



(A.9) 



with 



V2 \ ± 



± 



. The expansion of such a GEVP in e was written down in 
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with the intention that e is given by the HQET expansion parameter. We here jus ; 
use the solution. Its first order term in e yields the desired matrix element in the form 
eq. (3.1) in terms of the generalized eigenvectors Un of the lowest order (e = 0) GEVP 
of size X in a single channel A. 



See for example 21 , sect. 3.2. 
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